Fun Image and Signal Processing

Markus van Almsick

Wolfram Research

Fun Image and Signal Processing

Enhance

Due to Hollywood many people have a strange idea about what image processing can achieve.
YouTube “Enhance”  “Enhance even more”

An unrealistic example:

kernel=DiskMatrix[9]/Total[DiskMatrix[9],2];

img=ImageConvolve[,kernel]

ImageDeconvolve[img,kernel, Method->"SteepestDescent"]

Getting the point spread function wrong:

kernel2=DiskMatrix[10]/Total[DiskMatrix[10],2];

ImageDeconvolve[img,kernel2, Method->"SteepestDescent"]

Reality comes with noise:

img2=ImageEffect[img,{"GaussianNoise",1/128.}]

ImageDeconvolve[img2,kernel, Method->"SteepestDescent"]

Enhance

Find the Reflection

A typical family picture that needs some enhancement.

A closeup:

img=

To undo the bokeh, identify the point spread function given by the flash reflection in the eyes:

eyeCoordinates = {{376.5`,344.5`},{105.5`,362.5`}};

eyeReflections=Map[ImageTrim[img,{#},12]&,eyeCoordinates]

Taking the median of both point spread functions and all color channels:

eye1=First@eyeReflections;
Magnify[eye1,8]

eye2=RemoveAlphaChannel[ImageAlign@@eyeReflections];
Magnify[eye2,8]

kernel=ImageAdjust@ImageApply[Median,ColorCombine[{eye1,eye2}]];
Magnify[kernel,8]

Masking the point spread function by a blurred disk of radius 9:

kernel=ImageCrop[ImageMultiply[kernel,Blur@Image[DiskMatrix[9,ImageDimensions[kernel]],"Real"]],2 9+3];
Magnify[kernel,8]

Normalizing the kernel:

kernel=ImageMultiply[kernel,1/ImageMeasurements[kernel,"Total"]];

ListPlot3D[ImageData[kernel],MeshFunctions->{#3&},PlotRange->All]

Enhance

Deconvolve

deconvolved=ImageDeconvolve[img,kernel, Method->{"RichardsonLucy", "Preconditioned"-> False},Padding->"Reversed",MaxIterations->16]

Regularization

testResult=PeronaMalikFilter[deconvolved,6,0.01]

Comparison

TabView[{img,ImagePerspectiveTransformation[deconvolved,TranslationTransform[{1,1}],DataRange-> Full]}]

Fashion

Coloring Textile Patterns

Textile pattern:

pattern =

Pattern segmentation:

segmentation=ClusteringComponents[pattern,8];

Colorize[segmentation]

Fashion

Taking the Color Compositions from Nature

Quince

photo = ;

rules=MapIndexed[First[#2]->#1&,DominantColors[photo,Max[segmentation]]]

Row[{Style[photo,ImageSizeMultipliers-> {1,1}],"   ->   ", Style[Colorize[segmentation,ColorRules->rules],ImageSizeMultipliers-> {1,1}]}]

Fashion

Taking the Color Compositions from Nature

Autumn Forest

photo = ;

rules=MapIndexed[First[#2]->#1&,DominantColors[photo,Max[segmentation]]]

Row[{Style[photo,ImageSizeMultipliers-> {1,1}],"   ->   ", Style[Colorize[segmentation,ColorRules->rules],ImageSizeMultipliers-> {1,1}]}]

Fashion

Taking the Color Compositions from Nature

St. Petersburg

photo = ;

rules=MapIndexed[First[#2]->#1&,DominantColors[photo,Max[segmentation]]]

Row[{Style[photo,ImageSizeMultipliers-> {1,1}],"   ->   ", Style[Colorize[segmentation,ColorRules->rules],ImageSizeMultipliers-> {1,1}]}]

Fashion

Taking the Color Compositions from Nature

Chinese Door in the Forbidden City

photo = ;

rules=MapIndexed[First[#2]->#1&,DominantColors[photo,Max[segmentation]]]

Row[{Style[photo,ImageSizeMultipliers-> {1,1}],"   ->   ", Style[Colorize[segmentation,ColorRules->rules],ImageSizeMultipliers-> {1,1}]}]

Fashion

Taking the Color Compositions from Nature

North Sea Beach

photo = ;

rules=MapIndexed[First[#2]->#1&,RandomSample[DominantColors[photo,Max[segmentation]],Max[segmentation]]]

Row[{Style[photo,ImageSizeMultipliers-> {1,1}],"   ->   ", Style[Colorize[segmentation,ColorRules->rules],ImageSizeMultipliers-> {1,1}]}]

Fashion -> Art

Andy Warhol’s Marylin

img = ;

Separate color from grayscale for segmentation:

LabImages=ColorSeparate[img,"LAB"]

Create a color-only image:

MeanShiftFilter[ColorCombine[Join[{ConstantImage[0.5,{590,600}]},Take[LabImages,-2]],"LAB"],3,1/8]

Segement color distribution:

segmentation=ClusteringComponents[%,6];

Extract dominant colors from a grafiti:

rules=MapIndexed[First[#2]->#1&,DominantColors[,6]]

Rendering the recolored version:

ImageMultiply[Colorize[segmentation,ColorRules->rules],First[LabImages],2]

Speeding

Principles of the Doppler Effect

Sound sample of a Doppler effect

snd= ;

sampleRate=(QuantityUnit[8192]);

Speeding

Measuring the Frequency Shift

Apply a continuous wavelet transformation with a Gabor wavelet:

cwt=ContinuousWaveletTransform[snd,GaborWavelet[24],{6,24}];

WaveletScalogram[cwt,{2|3|4|5|6,_}]

Antialiased sub-sampling of the wavelet data to 32Hz:

subsamplingFactor = sampleRate/(QuantityUnit[32]);

waveletData2=
ArrayResample[
GaussianFilter[Abs@Map[First,cwt[All,"Values"]],Round[{{0,subsamplingFactor},{0,subsamplingFactor/4}}]],
{All,Scaled[1/subsamplingFactor]}
];

Measuring the scale shift:

scaleShifts=
Accumulate[Table[
FindPeaks[ListCorrelate[waveletData2[[All,t]],ArrayPad[waveletData2[[All,t+1]],4,0]], InterpolationOrder->3][[1,1]]-5,
{t,1,Last@Dimensions[waveletData2]-2}
]];

Converting shifts in wavelet scale back to change in frequency:

scaleToFrequency[cwd_]:=Evaluate[((cwd["SampleRate"])/(cwd["WaveletScale"])) (2^(-(#)/cwd["Voices"]))]&;

frequencyFactor=Normalize[scaleToFrequency[cwt][scaleShifts],Mean];

ListLinePlot[frequencyFactor]

Speeding

Fitting the Doppler Formula

Fitting Doppler model to measured change in frequency:

nlm=NonlinearModelFit[frequencyFactor,1-(((t/32-t0) (v^2))/(343.2(Sqrt[(d^2)+(((t/32-t0))^2) (v^2)]))),{A,v,d,t0},t]

nlm["ParameterTable"]

Result:
car speed = 19.5 m/s
min. distance = 4.7m  (at time 2.51693s)

UnitConvert[Quantity[nlm["ParameterTableEntries"][[2,1]],"m/s"],"mph"]

UnitConvert[Quantity[nlm["ParameterTableEntries"][[2,1]],"m/s"],"km/h"]

Blush

Heart Rate Estimation from Video

video=Image3DSlices[];

View the imported video:

ListAnimate[video]

Import frame acquisition times:

videoTs=0.001*{0,38,93,150,206,263,319,376,431,487,544,601,657,714,768,825,880,937,992,1049,1104,1160,1217,1273,1329,1385,1442,1498,1554,1611,1667,1723,1836,1893,1948,2004,2059,2116,2171,2227,2284,2340,2398,2453,2509,2567,2622,2678,2735,2791,2846,2903,2958,3015,3070,3127,3183,3239,3295,3351,3407,3464,3519,3574,3631,3686,3743,3799,3855,3911,3967,4024,4079,4135,4192,4248,4305,4416,4473,4529,4585,4641,4697,4753,4809,4864,4922,4977,5034,5089,5147,5203,5258,5314,5370,5427,5483,5540,5595,5651,5708,5764,5820,5877,5933,5989,6045,6102,6157,6214,6270,6327,6383,6438,6495,6551,6606,6662,6718,6831,6888,6943,7000,7055,7112,7167,7223,7280,7336,7393,7448,7504,7561,7616,7672,7729,7786,7843,7898,7955,8010,8066,8123,8180,8236,8292,8348,8403,8461,8517,8573,8629,8686,8741,8798,8855,8910,8967,9022,9079,9134,9191,9247,9303,9360,9417,9473,9585,9640,9696,9752,9808,9864,9922,9977,10034,10090,10145,10203,10258};

Blush

Face Detection

Determine the bounding box of the face:

faceBox=First@FindFaces[First@video]

HighlightImage[First@video,Graphics[{EdgeForm[Orange],FaceForm[],Rectangle@@faceBox}]]

faceBoxes=Map[First@FindFaces[#]&,video];

Find the median bounding box size:

avgFaceBoxSize=Round[{1/GoldenRatio,1} Map[Median,Transpose@Map[Differences,faceBoxes][[All,1]]],2]

Determine the regularized motion of the bounding box:

avgFaceBoxPositions=Transpose@Map[GaussianFilter[#,4]&,Transpose@Map[Mean,faceBoxes]];

Trim the video to the face bounding-box:

faces = MapThread[ImageTrim,{video,Map[{#-avgFaceBoxSize/2,#+avgFaceBoxSize/2}&,avgFaceBoxPositions]}];

Blush

Frame Registration

Determine the face-shifts with respect to the first frame:

faceShifts=Map[Last@FindGeometricTransform[First@faces,#,Method-> "Linear",TransformationClass-> "Translation"]&,faces];

Stabilize the video by undoing the shifts:

regFaces=MapThread[ImagePerspectiveTransformation[#1,#2, DataRange-> All,Resampling-> "Cubic"]&,{faces,faceShifts}];

Blush

Skin Detection

Define a typical skin classifier based on Lab color to detect facial skin:

skinClassifier= Compile[{{Lab,_Real,1}},Exp[-(((Lab[[2]]-0.25))^2)/(2 (0.06^2))] Exp[-(((Lab[[3]]-0.16))^2)/(2 (0.04^2))]];

Convert faces into Lab-color-space:

LabFaces = Map[ColorConvert[#,"LAB"]&,regFaces];

Determine the average skin probability distribution in the face bounding box:

skinWeight=ImageAdjust[ImageAdd@@Map[ImageApply[skinClassifier,#]&,LabFaces]]

Blush

De-Mixing Colors

For each frame, extract from all pixels the average Lab colors weighted by the above skin probability distribution:

LabSignal = Map[ImageMeasurements[#,"Mean",Masking-> skinWeight]&,LabFaces];

Taking a closer look at the color signals a and b:

ListLinePlot[Rest@Transpose@LabSignal,PlotLegends-> {"a","b"}]

Find the optimal de-mixing angle α by searching for a signal with the least differential volatility.

opt=Last@FindMinimum[Total[((Differences[{0,Cos[α],Sin[α]}.Transpose[LabSignal]])^2)],{α,0}]

Extracting the optimal pulse signal:

signal =({0,Cos[α],Sin[α]}/.opt).Transpose[LabSignal];

ListLinePlot[signal,Axes-> {Automatic,None}]

Blush

Pulse Measurement

Regularize the pulse signal with a bandpass filter that selects frequencies between 0.5Hz and 3Hz:

filteredSignal=BandpassFilter[signal,2π/18{0.5,3},11];

ListLinePlot[filteredSignal,Axes-> {Automatic,None}]

Extract the heart beats with FindPeaks:

heartBeats=FindPeaks[filteredSignal,2,InterpolationOrder-> 3, Padding-> 0.01]

ListLinePlot[filteredSignal,Epilog-> {Map[Line[{{First[#],0},#}]&,heartBeats]}, Axes-> {Automatic,None}]

Convert the frame numbers into times:

timeFunction=ListInterpolation[videoTs,InterpolationOrder-> 1];
Ts=Map[timeFunction,heartBeats[[All,1]]]

Extract heart beat intervals:

ΔTs=Differences[Ts]

Median heart rate and its deviation:

mΔT=Quantity[Median[ΔTs],"Seconds"]

dΔT=Quantity[MedianDeviation[ΔTs],"Seconds"]

Median pulse:

Quantity[60,"Seconds"]/mΔT

The End

Any Questions?